AD-A280  990 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


V5  ^ 


THESIS 


-rir* 

1 1\,-' 


^-ItELECTEf'l 

juLaami 


AN  INVESTIGATION  OF 

THE  TRANSONIC  PRESSURE  DRAG  COEFFICIENT 
FOR  AXI-SYMMETRIC  BODIES 

by 

Eddy  Priyono 
March,  1994 

Thesis  Advisor;  Oscar  Biblarz 

Approved  for  public  release;  distribution  is  unlimited. 


94-20400 


94  7  5  113 


DISCLAIMEI  NOTICE 


TfflS  DOCUMENT  IS  BEST 
QUALITY  AVAILABLE.  THE  COPY 
FURNISHED  TO  DTIC  CONTAINED 
A  SIGNffICANT  NUMBER  OF 
COLOR  PAGES  WHICH  DO  NOT 
REPRODUCE  LEGIBLY  ON  BLACK 
AND  WHITE  MICROFICHE. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  0MB  No.  0704 


Public  rapoiting  burdeo  for  thii  coUection  of  iafoniiatioa  ii  MliiiiaMd  to  avnafo  1  hour  per  leipoiue,  including  the  time  for  reviewing  iniuuction, 
searching  existing  dau  souroea,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  coUection  of  information.  Send  commenu 
regarding  this  burden  estimate  or  any  other  aspect  of  this  coUection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washingon 
headquarters  Services,  Directorale  for  Information  Operations  and  Raportt,  1213  Jefferson  Davis  Highway,  Suite  1204.  Arlington,  VA  22202-4302,  and 
to  the  OfRce  of  Manageoiem  and  Budget,  Paperwork  Reduction  Projea  (0704-0188)  Washington  DC  20303. 

1.  AGENCY  USE  ONLY  (Leave  blanki  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

March  1994.  Master’s  Thesis 

4.  TITLE  AND  SUBTITLE  AN  INVESTIGATION  OF  THE 

TRANSONIC  PRESSURE  DRAG  COEFFICIENT  FOR  AXI- 
SYMMETRIC  BODIES 

5.  FUNDING  NUMBERS 

6.  AUTHOR(S)  Eddy  Priyono 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Postgraduate  School 

Monterey  CA  93943-5000 

8.  PERFORMING 
ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10. 

SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not 
reflect  the  official  policy  or  position  of  the  Department  of  Defense  or  the  U.S.  Government. 

12a.  DISTRIBUTION/ AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

12b.  DISTRIBUTION  CODE 
•A 

13. 

ABSTRACT  (meudmum  200  words) 

This  thesis  investigates  the  pressure  drag  coefficient  in  the  transonic  regime  over  an  axi-synunetric  body,  with  a 
set  of  unique  contour  surfaces  developed  in  a  previous  thesis.  The  contour  surfaces  were  obtained  by  an  exact 
solution  of  the  small  perturbation  transonic  equation,  using  the  guidelines  and  tools  developed  at  NFS.  In  this  work, 
Computational  Fluid  Dynamics  (CFD)  was  not  only  used  to  compute  the  afterbody  contour  surface,  but  also  to 
investigate  a  conical  afterbody  and  complete  bodies,  which  are  composed  of  an  arbitrary  forebody  (ellipsoid)  and  a 
variable  afterbody  (contour  and  conical).  Euler  as  well  as  Navier-Stokes  flow-solvers  were  applied  to  the  geometries 
of  interest,  giving  Mach-number  contours  for  viscous  and  inviscid  flow,  pressure  drag  coefficient  magnitude,  and 
depicting  shock  wave  location.  On  the  basis  of  these  results,  it  can  be  verified  that  our  contour  surface  afterbodies 
will  decrease  by  15%  the  peak  of  the  pressure  drag  coefficient  (C,^  versus  Mach  number  curves  in  the  transonic 
regime.  These  results  can  be  used  to  design  low  pressure  drag  surfaces  for  such  as  missiles,  projectiles  and  aircraft 
engine  nacelles. 

14.  SUBJECT  TERMS  Tnuisonic  pressure  drag  coefficient  for  axi-symmetric 
bodies. 

15. 

NUMBER  OF 
PAGES  94 

16. 

PRICE  CODE 

17. 

SECURITY  CLASSin- 
CATION  OF  REPORT 
Unclassified 

18. 

SECURITY  CLASSin- 
CATION  OF  THIS  PAGE 
Unclassified 

19. 

SECURITY  CLASSin- 
CATION  OF 

ABSTRACT 

Unclassified 

20. 

UMITATION  OF 
ABSTRACT 

UL 

NSN  7540-01-280-5500  Standard  t^orin  298  (Rev.  2-89) 


Prescribed  by  ANSI  Std.  239-18 


1 


Author: 


Approved  by: 


Approved  for  public  release;  distribution  is  unlimited. 

AN  INVESTIGATION  OF 

THE  TRANSONIC  PRESSURE  DRAG  COEFFICIENT 
FOR  AXI-SYMMETRIC  BODIES 

by 

Eddy  Priyono 

Major,  Indonesian  Air  Force 
Mechanical  Engineer, 

Institute  Technology  of  Surabaya,  1979 

Submitted  in  partial  fulfillment 
of  the  requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  AERONAUTICAL  ENGINEERING 


from  the 


Department  of  Aeronautics  and  Astronautics 


II 


ABSTRACT 


This  thesis  investigates  the  pressure  drag  coefficient  in 
the  transonic  regime  over  an  cuci- symmetric  body,  with  a  set  of 
unique  contour  surfaces  developed  in  a  previous  thesis.  The 
contour  surfaces  were  obtained  by  an  exact  solution  of  the 
small  perturbation  transonic  equation,  using  the  guidelines 
and  tools  developed  at  NPS.  In  this  work,  Computational  Fluid 
Dynamics  (CFD)  was  not  only  used  to  compute  the  afterbody 
contour  surface,  but  also  to  investigate  a  conical  afterbody 
and  complete  bodies,  which  are  con^osed  of  an  arbitrary 
forebody  (ellipsoid)  and  variable  afterbody  (contour  and 
conical) .  Euler  as  well  as  Navier- Stokes  flow- solvers  were 
applied  to  the  geometries  of  interest,  giving  Mach-number 
contours  for  viscous  and  inviscid  flow,  pressure  drag 
coefficient  magnitude,  and  depicting  shock  wave  locations.  On 
the  basis  of  these  results,  it  can  be  verified  that  our 
contour  surface  afterbodies  will  decrease  by  15%  the  peak  of 
the  pressure  drag  coefficient  (C^)  versus  Mach  number  curves 
in  the  transonic  regime.  These  results  can  be  used  to  design 
low  pressure  drag  surfaces  for  such  shapes  as  missiles, 
projectiles  and  aircraft  engine  nacelles. 
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Z.  ZHTRODUCTZOH 


The  flight  of  an  object  over  a  wide  range  of  speeds  has  a 
critical  transition  zone  where  both  subsonic  and  supersonic 
types  of  flow  exist.  This  speed  regime  is  referred  to  as  the 
transonic  range.  The  critical  aerodynamic  behavior  occurs  in 
the  range  0.8  <  M  <  1.2  depending  on  the  object,  where  the 
aerodynamic  coefficients  have  been  found  to  change  by  as  much 
as  100%.  The  behavior  of  the  aerodynamic  force  components  is 
usually  characterized  by  a  rapid  increase  in  the  coefficients 
followed  by  a  sharp  drop,  in  other  words,  a  peak  value  for  the 
coefficients  arises  in  this  regime  [Ref.l  and  2]. 

One  of  the  aerodynamic  components,  aerodynamic  drag, 
represents  a  significant  adverse  force  on  all  flying  objects 
such  as  aircraft,  missiles  and  projectiles.  A  high  drag  force 
reduces  the  craft's  range  capability  or  equivalently  requires 
more  energy  to  achieve  a  certain  range.  Any  effort  to  reduce 
the  drag  coefficient  in  the  design  process  must  concentrate  on 
reducing  the  wake  and  pressure  drag  (inclusive  of  wave  drag) 
contributions  to  the  total  drag. 

A  traditional  approach  to  investigate  aerodynamic 
characteristics  is  based  on  wind  tunnel  data  and  actual  flight 
testing  as  a  source  for  improvement  of  configurations  to  get 
the  best  results.  Unfortunately,  both  wind  tunnel  and  flight 
tests  are  considerably  expensive  and  time  consuming. 
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particularly  in  the  transonic  flow  regime.  Because  of  the 
great  lateral  influence  of  profiles  at  transonic  speeds  (shock 
reflection  of  the  wall  bounces  back  to  the  model)  ,  models  must 
be  made  extremely  small  compared  with  wind  tunnel  dimensions 
and  this  introduces  great  experimental  difficulties. 
Consequently,  the  data  obtained  at  transonic  speeds  is 
considerably  less  reliable  than  at  either  subsonic  or 
supersonic  speeds.  In  addition,  from  a  mathematical  point  of 
view,  even  the  two  dimensional  small  perturbation  potential 
equation  for  transonic  flow  retains  one  non-linear  term  which 
is  essential  for  non-divergent  solutions  at  Mach  one,  but  this 
non-linear  partial  differential  equation  has  proven  to  be 
difficult  to  solve.  Such  inherent  difficulties,  coupled  with 
the  presence  of  shocks  in  the  flow  which  cause  boundary  layer 
separation,  have  resulted  in  the  creation  of  many  approximate 
methods  of  solution  which  are  employed  in  the  design  of 
transonic  airfoils  and  the  like  [Ref. 3,  4  and  5]  .  On  the  other 
hand,  the  use  of  niomerical  simulation  known  as  Computational 
Fluid  Dyncunics  (CFD)  to  predict  aerodynamic  characteristics 
greatly  increases  possibilities  to  improve  design  optimization 
at  relatively  low  cost  and  allows  for  ease  of  design  changes. 

Finally,  using  the  latest  capabilities  of  Euler  as  well  as 
Navier- Stokes  flow- solvers,  it  has  been  possible  to  compute 
the  flow  over  axi- symmetric  bodies  with  various  contours  in 
the  transonic  regime. 
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II  TRANSONIC  FLOW 


Transonic  flows  are  characterized  by  the  simultaneous 
presence  within  the  flow  field  of  both  subsonic  and  supersonic 
regions.  The  properties  of  transonic  flows  can  be  derived  from 
the  general  equations  of  gas  dynamics,  namely,  the  equations 
of  state,  continuity,  momentum  and  energy.  The  following 
derivation  based  on  assumptions  that  the  flow  is  steady, 
irrotational  and  isentropic  with  no  energy  transfer,  no  body 
forces  and  no  shear  stresses  (inviscid  flow) . 


1.  Small  Perturbation  Theory 

Starting  with  the  equation  of  motion  for  steady, 
isentropic,  inviscid  flow  in  the  index  notation  form  [Ref  .4]  : 


du.  2  du^. 

^dxj  dxj, 


(1) 


Small -Perturbation  Theory  gives  the  velocity  field  as 

=  U  +  u  ;  Uj  =  V  ;  U3  =  w  (2) 

where  U  denotes  the  free  stream  velocity  in  the  x  direction 
and  u,  V,  w  are  called  perturbation  velocity  components  in  the 
direction  x,  y,  z  respectively. 
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Writing  eq. (1)  out  in  full,  and  substituting  the  velocity 
field  eq. (2) ,  gives  the  equation  in  terms  of  perturbation 
velocities. 


\  dx  ay  dzf  dx  dy  oz 


dv  .  dw 


du 


dv 


dw 


(3) 


*  t) '  if  *  I?)  * *  i) 


dv  ^  dw\ 


dw  ^  du\ 


From  the  energy  equation  for  a  perfect  gas,  the  speed  of  sound 
(a)  can  be  expressed  in  terms  of  the  perturbation  velocities. 


(U+u)^  +  V*  +  ^  ^ 

2  *  (Y-l)  "  T  ^  (Y-l) 


(4) 


or 


=  al 


_  (y-1) 
2 


{2uU 


u‘ 


w^) 


(5) 


Substituting  eq. (5)  into  eq. (3) ,  dividing  by  a„^,  and 
rearranging  the  terms,  gives  the  full  exact  equation  in  terms 
of  perturbation  velocities  and  free  stream  Mach  number.  This 
equation  contains  linear  terms  on  the  left-hand  side,  but  on 
the  right-hand  side  the  terms  are  nonlinear. 
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n-M^) 


du 

dx 


—  +  ^ 
dy  dz 


(Y+1) 


u 

U 


+ 


+ 


(Y+1) 

(Y-1) 

( V-2+nr2) 

du 

2 

2 

dx 

(Y+1) 

(Y-1) 

(  v2  +  u2) 

dv 

2 

2 

dy 

(Y+1) 

+ 

(Y-1) 

iu^+v^) 

dw 

2 

2 

dz 

(61 


vwl  dw^  dv 
c;2\3y  dz 


If  the  perturbation  velocities  are  small  (u/U,v/U,w/U  <<  1) , 
eq. (6)  can  be  simplified  by  neglecting  the  terms  containing 
squares  of  the  perturbation  velocities  on  the  right  hand  side, 
yielding  ; 
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For  further  simplification,  in  eq. (7)  all  the  terms  on  the 
right  hand  side  can  be  neglected,  in  comparison  to  those  on 
the  left  hand  side.  This  gives  the  linear  equation,  which 
contains  only  perturbation  velocities  and  is  valid  only  for 
subsonic  and  supersonic  flow. 


dx 


dv 

dy 


(8) 


For  transonic  flow,  where  1,  the  coefficient  of  du/dx  on 

the  left  hand  side  becomes  very  small,  but  it  is  not  correct 
CO  neglect  the  first  term  on  the  right  hand  side  of  eq. (7) . 
Therefore,  the  governing  equation  of  transonic  f3ow  in  term  of 
perturbation  velocities  is  as  follows: 


(1-^d) 


du 


dv  ^  dw 

~Uy  ~3z 


wf (7+1) 


u  du 
U~^ 


(9) 


For  irrotational  flow,  a  perturbation  velocity  potential 
exists. 


u 


(10) 


Substitution  of  eq. (10)  into  eq. (9) ,  gives  the  governing 
equation  for  transonic  flow  in  terms  of  the  velocity 
potential . 


6 


(11) 


‘'■'^’ aP  IP  IP  —SJijp 

For  bodies  of  revolution,  it  is  convenient  to  use  cylindrical 
coordinates  {x,r,d)  where  x  is  aligned  with  the  body  axis.  The 
velocity  components  corresponding  to  (x,r,ff)  are  Uj^,  u^.  and 
U0,  respectively.  The  velocity  potential  in  the  cylindrica 
coordinates , 

Transforming  Cartesian  coordinate  eq. (9)  into  cylindrical 
coordinate,  gives  the  governing  equation  of  transonic  flow. 

d^<t>  ^  d^4>  ^  1  ^  1  d^<t>  =  d4>  d^4>  (13) 

‘  3P  IP  PIP 

For  axially  symmetric  flow,  where  the  conditions  are  the  same 
in  every  meridian  plane,  there  is  no  variation  with  d,  so  the 
small  perturbation,  non-linear,  auci- symmetric  transonic 
potential  equation  can  be  written  as  follows  : 

d^(t>  ^  (14) 

dr2  rTr  9^2  u  ~Sx  dx^ 
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Or  rewriting  in  shorthand  notation, 


<t> 


rr 


XX 


(15) 


Reference  6  introduces  a  modified  potential  equation  for  eoci- 
symmetric  flows 

where  the  modified  velocity  potentials  are  : 

♦jc  -  (T-l)  ^  (17) 


2.  An  exact  solution  for  axi- symmetric,  transonic  flow. 

Solutions  to  the  modified  transonic  eq. (16)  have  been 
given  by  Biblarz  [Ref. 6  and  7]  by  using  the  separation  of 
variables  approach  with  a  potential  function  <^(x,r)  of  the 
form 

Substituting  the  above  <t>  function  in  the  modified  transonic 
eq. (16)  results  in  two  ordinary,  second  order,  non-linear 
differential  equations 
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(19) 


dj 

dx^ 


-  =  0 


and 


dr^ 


1  dll 

TgE 


-  Xij2  _  0 


where  X  is  the  separation  constemt. 

The  solution  to  the  first  differential  eq. (19)  is 
by  multiplying  both  sides  df/dx. 


or 


Thus 


dx 


'3x  fjx2 


-§(X{)  =0 


*  0 


+  a 


1 


where  a  is  a  constant  given  by  Ref. 6. 
Rearranging  eg. (23) 


(20) 


obtained 


(21) 


(22) 


(23) 


(24) 
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and  integrating  eq. (24)  gives. 


(25) 


where  a  and  Xq  are  integration  constants. 

The  solution  to  the  second  non-linear  ordinary 
differential  equation  (20) ,  is  obtained  by  an  outer  expansion 
method.  [Ref. 8] 

fl(r)  =  -i-  (l-^^)fl(r)  *  (l-^^)2f2(r)  +  ...  (26) 

Where  (l-Mn^)  represents  a  small  parameter  and  the  first  term 
is  the  purely  sonic  solution. 

By  taking  first  and  second  derivatives  eq.  (26)  cuid 
subtituting  them  into  eq. (20)  yields, 


n  (r)  * 


1 

Xr2 


1-/1^  laXr^'^'^2) 


(1-Mf  )^^2\2y2(yF-»2)  ^ 

28+8>/i’ 


k27) 


Al-Hashel  [Ref. 9]  reported  on  eq. (26)  and  eq. (27) , 
inqplemented  the  boundary  conditions  with  the  constants  a,  C^, 
a  and  X,  gave  the  final  results  in  the  new  varied)les, 


f 


3X 

TUI 


1 

■J 


(28) 
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^wer*2) 


r(yF-2)  X  a 


(29) 


X  m 


(30) 


V 


n  X 

(a 


(31) 


then 


Si 


_ d| _ 

[p*  u-idi'-’”*]’ 


(32) 


where  Xq-0 ,  and 


V 


50.63 


^2WB-*2) 


(33) 


Equation  (32)  has  been  numerically  integrated  and  plotted  on 
figure  (l)  as  f  (x)  versus  x,  and  eq.  (33)  will  be  evaluated  and 
plotted  in  figure  (2)  as  ^(f)  versus  f  for  =  1.05,  1.1  and 
1.2.  A  "patching"  technique  discussed  in  Ref. 9  has  been  used 
here. 
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III.  PRESSURE  COEFFICIEET  AND  BOUNDARY  SURFACES 


1.  Pressure  Coefficient 

Liepmann  and  Roshko  [Ref  .4]  define  the  pressure  coefficient 
as, 


=  =  2 
0.5p«Ej2  yM^ 


(34) 


From  the  isentropic  relation,  we  have  pressure  ratio  in  terms 
of  Mach  number  and  after  substituting  into  eq. (34)  yields 


2  +  {7-l)i>^ 


17/ W-i) 


(35) 


Introducing  m;^,^  =  =  u^/a^  and  using  energy  equation 
(recall  eq.4  and  eq.5),  the  pressure  coefficient  can  be 
expressed  in  terms  of  the  perturbation  velocities. 
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Using  the  binomial  expemsion  on  the  expression  inside  the 
square  brackets  in  eq.  (36)  ,  we  obtain  the  pressure  coefficient 
in  the  form 


IR. 

u 


(37) 


For  axi- symmetric  flow,  in  cylindrical  coordinates  where  u=Ujj 
and  (v^+w^)  =  u^.^  ,  substituting  into  eq.  (37)  yields 


-*•  (l-Ad)f-J 

L)  (i 

1 

(38) 


The  linearized  pressure  coefficient  approximation  for  suci- 
symmetric  flow  turns  out  to  be 


(39) 


Recall  the  modified  axial  velocity  potential,  eq. (17) , 


<t>x  = 


(40) 


thus 


^ 

Mi(7+1) 


(41) 
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Sxibstitute  eq.  (41)  into  eq.  (39)  ,  yields 


•  = 

^  «^(7-^l) 


(42) 


The  derivative  of  the  potential  function  eq. (18)  with  respect 
to  X  becomes 

=  n  ^  ^  n-ti)  (43) 

Ref. 6  introduces  relation  of  the  constants  Cj^,  a,  a  and  X  in 
expression. 


a  =  ±Cill-ii  |^•■^574 


(44) 


and 


Cl]  ^1.2*26 
"XT  ;^0.7574 


1.08x10"2 


(45) 


Rewriting  eq. (23)  and  substitute  the  constant  a  eq. (44) 
becomes 


di 

■asE 


3X>2 
2  ^ 


Cii  ii-ii)  I 


1.7574 


1 

1 


(46) 
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Recall  the  new  variable  eq. (28) 


f  =  I 


3X  1-5 


Substituting  eq.  (47)  into  eq.  (46)  and  factoring  Cj^  out,  yields 


4  = 


Recall  and  arrange  new  varicUale  eq.  (31) 


«  .  (a  X) 0-4142 
V  -  - c -  V 


Rewriting  eq. (43)  the  modified  potential  function  and 
inserting  eq. (48)  and  eq. (49) ,  becomes 


=  (aX)°;4^42^  ^3  I  1.7574]?  ^  (1-wf) 


aO -4142^7 


f  m  — 


(i-wf) 


15 


Substituting  the  expression  of  constants  eq. (45)  into  eq. (51) 


= 


0.2208  ^  V  *  \  (l-^^)  |^•'^574 


1 

T 


(52) 


Finally,  rewriting  the  pressure  coefficient  Cp  as 


= 


-2 


—  b.2208^fp..|  (l-^^) 

•  1)  L  L  j 


(53) 


2 .  Boimdary  Surfaces 

For  an  inviscid  flow,  the  condition  to  be  applied  at  the 
surface  of  a  solid  boundary  is  that  the  direction  of  the  flow 
velocity  vector  must  be  tangent  to  the  solid  surface.  In  other 
words,  the  velocity  vector  is  everywhere  at  right  angles  to 
the  normal  to  the  solid  boundary  [Ref. 10].  In  addition  the 
boundary  condition  requires  that  the  gradient  of  potential  0 
vanish  far  ahead  of  the  body.  In  terms  of  perturbation 
velocities  this  boundary  condition  becomes 


l^\  =  Jfi 

\  djc/  sur/«C» 


(54) 
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Recall  the  modified  velocity  perturbation  potential  eq. (17) 


^OD 


(55) 


thus 


(56) 


substituting  eq. (56)  into  eq. (54) , 


/  dr\  _ 

(7+1) 


(57) 


Taking  the  d.*  fferential  of  eq.  (18)  with  respect  to  r  gives 


-  t  dv 
^  dr 


(58) 


Expressing  eq. (58)  in  terms  of  new  variables,  recall  eq. (28) , 
eq. (29)  and  eq. (31) 

<t>j,  =  0.0848  I  ^  (59) 
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Recalling  eq.  (33)  euid  taJcing  the  derivative  with  respect  to  f  , 
then  substituting  eq. (59)  into  eq. (57) ,  yields 


(^\  —1+2.8284  “M+0.1512  (1-Wf)2  f®  *” 


(60) 


Replacing  the  left  hand  aide  of  eq.  (60)  with  the  new  variables 
eq. (29)  and  eq(30),  then  rearranging  it  becomes 


^  =  _0-0326  ^  -I£+2.8284  .  1512  (1-W(^)  2  f6.657 


(61) 


For  further  arrangement,  by  separating  variables  eq. (61) , 
yields 


df 


-p-2.828A\l-M!i\  i'i-«2««-0.1512  (l-/yf)2  rs-ss? 

I  dSt 


(62) 


Then,  integrating  both  sides  of  eq. (62)  becomes 
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* 

/ 


df 


-?--2.8284|l-/d;2|f  -0.1512  (1-wf)^  f*  ”’ 


-0.0326 


(63) 


where 


1.2074 

l-^ff  |‘>-2071 


(64) 


Al-hashel  [Ref. 101  has  developed  and  computed  eq. (63) 
using  numerical  integration,  to  determine  the  boundary 
surfaces  in  dimensional  and  non-dimensional  (normalized)  fom 
for  Mlj,  =  1.05,  1.10  and  1.20  as  depicted  in  Figure  (3)  and 
Figure  (4)  .  Based  on  these  calculations,  the  geometric  grid  of 
the  afterbody  for  Mach  number  1.10  and  1.20  are  developed  for 
further  study  with  Computational  Fluid  Dynamic  (CFD) . 

This  thesis  research  also  examines  conical  afterbodies  as 
a  solid  afterbody  boundary  surface  with  base  diameter  ratio 
(djj/d„)  of  0.50  and  0.75  and  conical  turning  angle  (jS)  of  26.6 
and  14.0  degree  respectively  [Ref  .11  and  12]  .  Then,  for 
further  investigation,  the  coii^lete  bodies  as  a  solid  boundary 
surface  are  generated,  with  a  kind  of  forebody  (ellipsoid) 
joined  with  the  contour  surface  afterbodies  a,  well  as  the 
conical  afterbodies. 
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IV.  DRAG 


Drag  is  one  of  the  aerodynamic  force  components  parallel 
to  free  stream  velocity  (U») .  It  represents  a  significant 
adverse  force  on  all  flying  objects.  Basically,  the  drag  force 
is  divided  into  two  categories,  the  drag  caused  by  forces 
acting  normal  to  the  boundary  surface  which  is  called  pressure 
drag  (inclusive  forebody,  base  and  wave  drag)  and  that  arising 
from  the  tangential  forces  acting  on  the  surface,  by  virtue  of 
viscosity  which  is  called  viscous  drag  or  skin  friction. 
Figure  (5)  shows  the  components  of  the  drag  as  function  of 
Mach  number  and  the  methods  used  to  compute  body  drag  in  four 
Mach  number  regimes  (Ref.SJ. 


Figure  5.  Components  of  the  Drag  [Ref. 5] 
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1 .  Skin  Friction 


Skin  friction  is  the  result  of  shearing  force  within  the 
boundary  layer  of  a  viscous  flow,  acting  tamgentially  to  a 
surface  in  motion  relative  to  the  fluid.  The  amount  of  viscous 
resistance  depends  upon  whether  the  flow  is  laminar  or 
turbulent  [Ref  .5]  .  Krasnov  [Ref  .13]  introduces  the  leuninar  euid 
turbulent  skin  friction  for  flat  plates,  based  on  the  boundary 
layer  theory,  for  con^ressible  flow.  These  formulas  also  valid 
for  bodies  of  revolution  with  infinite  length. 

Lcuninar  flow.  Re  <  10®  : 


^tL  ~ 


1.32 

(i?e)  2 


.  1 

(1  +  0.03/d) 


(65) 


Turbulent  flow.  Re  >  10®  ; 


Cgj,  = 


[logioi?©  (l+0.2A«f)-^-''«l'^-®® 

(1+0. 1«^)  ^ 


(66) 


where  Re  is  the  Reynold  number 
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2 .  Prttasure  Drag 

Liepmann  and  Roshko  [Ref. 4]  introduce  the  pressure  drag 
formula  (inclusive  wave  drag  and  base  drag)  for  eoci- symmetric 
bodies  as, 


L 

D  =  j  p  dS  -  PgS(L) 

0 

(67) 

L 

=  j  (p-pJdS  +  (p,-Pa)5(L) 

0 


or,  in  dimensionless  form. 


Cr, 


D 

q^(L) 


1 

S{L) 


/ 


C  ^dx 
p  dx 


+  C. 


pB 


(68) 


where  S{x)  =  irr^  is  the  cross-sectional  area  of  the  body  at  x. 

dS(x)  =  2nrdr  =  2xr-^dx  (69) 

dx 


or 


0^ 

dx 


=  2itr 


dx 


(70) 
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The  first  term  in  sc[uare  bracket  eq.  (68)  com  be  solved  by 
inserting  eq. (53) ,  eq. (60)  and  eq. (70)  yielding. 


-2 


id  (7+1) 


l).2208nlP  +  l  (1-id)  ‘^  +  (l-ld) 


(71) 


(2Trr) 


0.0326  I 


-^5 - «fl|+2.8284|l-ldl^^*®^®*+0.1512  (1-id) 


Wx 


Recalling  new  voiriables  eq.  (28) ,  eq.  (29),  eq.  (30)  amd 
inserting  0^  -  1-H.^  into  eq.  (71)  yields. 


I  ■ 


tt 

■ 

1 

X 

[ 

-2 

1 

p.2208ij(f+|/3®*®^*®|]  '^+j82 

» 

5.1302cJ 

i 

(l-j52)  (Y+l) 

j^0,7071  ^0.2071 

»  4 

(72) 


f 


0.0326 


(1-|82)  (y+1) 


I 


+2 . 8284 1  I  •  ®2®*  +0 . 1512j8*f  ® 


where  Cj^,  a  amd  X  are  constamts. 
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The  second  term  in  square  bracket  eq. (68)  is  the  contribution 
to  the  drag  by  the  base  pressure  pg.  Ref. 4  claims,  that  the 
values  of  the  base  pressure  coefficient  Cpg  must  be  obtained 
experimentally.  However  Krasnov  [Ref. 13]  introduces  the 
boattail  drag  coefficient  for  a  conical  afterbody  as. 


-pB 


=  0.002  k).8  +  ^ 


(73) 


where  : 

6  «  turning  angle  of  the  conical  body 

Sg  -  ratio  of  the  base  area  to  the  mid- section  area 

In  this  work,  the  skin  friction  coefficient  as  well  as  the 
pressure  drag  coefficient  of  the  geometries  of  interest,  at 
specified  Mach  number,  can  be  obtained  directly  from  the 
results  of  CFD  calculations.  The  Euler  solution  will  give  the 
pressure  drag  and  the  Navier- Stokes  solution  will  give  both. 
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V.  CGMPUTJLTZOKAL  FLTJZO  DZlIMflCS 


The  rapid  advancement  in  the  speed  of  con^uters  and  their 
enormous  memory  size  has  led  to  the  emergence  of  the  field  of 
computational  fluid  dynamics  (CFD) .  This  branch  of  fluid 
dynamics  conplements  its  experimental  emd  theoretical 
branches,  by  providing  an  alternate,  cost-effective  means  of 
simulating  real  flow.  It  also  offers  the  only  means  of 
excunining  theoretical  advances  for  conditions  unavailable 
experimentally.  As  a  model  based  method,  CFD  can  provide  the 
convenience  of  being  cQsle  to  switch  off  specific  terms  in  the 
governing  equations  [Ref. 14],  so  as  to  assist  the  researcher 
in  understanding  the  contributions  of  various  physical 
factors . 

In  this  work,  CFD  was  used  to  confute  the  cuci- symmetric 
flow  over  the  afterbody  geometry  of  models  only  (the  boundary 
surfaces  obtained  by  the  small  perturbation  method  [Ref  .9]  and 
conical  afterbody)  and  over  coo^lete  body  models  which  are 
composed  of  forebody  (ellipsoid)  and  afterbody. 

1 .  Grid  Generation 

The  conputer  programs  GRAPE  [Ref.  15]  euid  GRIDGEN2D 
[Ref  .16]  are  tools  used  to  generate  two-dimensional  structured 
grids  aisout  airfoils  cuid  other  shapes  by  the  use  of  algebraic 
or  Poisson  differential  equation  solvers.  GRAPE  was  used  for 
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the  geometry  with  smooth  contour  surfaces,  while  GRIDGEN2D  was 
used  for  the  geometry  with  a  conical  afterbody  with  non- smooth 
points.  Outer  and  inner  boundaries  were  specified  as  the  C- 
type  grid  for  afterbody  models  only,  while  the  0  type  grid  was 
for  the  complete  body  models,  where  both  type  of  grids  treat 
the  surface  of  the  body  as  the  inner  boundary.  The  important 
characteristics  in  a  grid  generation  technique  are  the  cibility 
to  specify  the  spacing  between  mesh  points  at  the  boundary,  in 
the  direction  normal  to  the  boxindary,  and  the  control  of  the 
angles  with  which  mesh  lines  intersect  the  boundaries  which  is 
known  as  orthogonality  [Ref .15] . 

Figure  (6)  is  a  typical  output  of  program  GRAPE  for  an 
afterbody  with  small  perturbation  solution  contour  for  Mach 
1.10  {SPS_l.l),  C-type  grid  with  the  grid  size  of  115x60,  and 
figure  (7)  is  the  complete  body  with  the  grid  size  of  152x60 
(0-type  grid) . 

To  develop  a  three-dimensional  grid  from  two-dimensional 
grid  (output  program  GRAPE),  the  FORTRAN  code  called  HALF.F 
(Appendix)  was  used  to  write  out  the  half  of  the  2D  grid  and 
convert  them  into  3D  plane -grid.  Then,  the  FORTRAN  code 
ROTATEGR.F  (J^pendix)  was  used  to  rotate  the  3D  plane -grid  for 
11  planes  plus  2  more  to  generate  a  3D  volume  grid  of  half  of 
an  cixi- symmetric  body  surface  as  shown  in  figure  (8)  and 
figure  (9)  with  grid  size  of  50x13x60  and  77x13x60 
respectively.  These  grids  are  ready  for  further  processing 
with  a  flow- solver. 
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Prograun  GRIDGEN2D  gives  the  pattern  of  how  to  build  2D- 
grid.  It  consists  of  four  subfaces  which  are  treated  as 
boundaries.  Grid  spacing  is  determined  by  setting  up  the 
distribution  points  of  each  pair  subfaces.  It  uses  equal 
spacing  along  the  body  surface,  while  the  direction  normal  to 
the  body  surface  uses  geometric  spacing  with  specified  width 
in  the  beginning.  Figure  (10)  and  Figure  (11)  are  typical 
outputs  of  GRIDGEN2D  for  an  afterbody  and  for  a  complete  body 
with  the  grid  size  of  58x60  and  77x60  respectively.  FORTRAN 
code  D2D3.F  (Appendix)  is  used  to  convert  2D  grid  output  from 
GRIDGEN2D  into  3D  plane -grid.  Then  we  apply  the  FORTRAN  code 
ROTATEGR.F  to  generate  a  3D  volume  grid  of  half  of  an  axi- 
symmetric  body  as  depicted  in  Figure  (12)  and  Figure  (13)  with 
the  grid  size  of  58x13x60  and  77x13x60. 

In  this  research,  we  also  attempted  to  develop  a  fine  grid 
for  complete  bodies  (0-type  grid) ,  where  the  radius  of  outer 
boundary  is  set  up  to  be  five  times  of  the  body's  length.  The 
steps  to  obtain  3D  volume  grid  are  similar  to  those  described 
above  with  the  final  grid  size  of  77x9x120  and  the  final 
geometry  of  a  quarter  of  an  axi- symmetric  body  instead  of  a 
half  body.  Figures  (14)  and  (15)  are  typical  fine  grids  of  the 
complete  body  with  a  conical  afterbody  and  a  complete  body 
with  contour  surface  afterbody. 
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2.  Flow- solver 


The  OVERFLOW  progrcun  [Ref.  17]  was  developed  by  NASA  Ames 
Research  Center.  It  uses  either  3-D  Euler  or  Navier-Stokes 
flow- solvers  for  inviscid/viscous  flow,  by  setting  True  or 
False  the  parameter  VISINP  (viscosity  input)  in  the  input  file 
(overflow. in) .  Before  flow-solver  code  (OVERFLOW)  is  applied 
on  the  grid  file,  a  formatted  3D  grid  file  named  "grid. for" 
must  be  converted  into  an  unformatted  grid  file  nauned 
"grid. in",  by  using  FORTRAN  code  called  READX.F  (Appendix). 
Then,  the  NAMELIST  input  file  parameter  specification  must  be 
written  for  running  OVERFLOW,  it  is  called  "overflow. in" .  The 
input  parameter  consists  of  the  number  of  iterations, 
timesteps,  calculation  methods,  smoothing,  type  of  flow  and 
boundary  conditions  for  each  grid.  The  value  of  angle  of 
attack  (ALPHA)  depends  on  the  orientation  of  the  grid  in  the 
coordinate  system.  In  this  case,  ALPHA  is  180°  (flow  comes 
from  the  x-positive  to  the  x-negative  direction)  .  The  boundary 
conditions  depend  on  the  geometry  corresponding  to  the  final 
3D  volume  grid.  Both  input  parameters,  the  angle  of  attack  and 
the  boundary  conditions  for  each  geometric  shape,  are 
tcd^ulated  in  Tcdjle  1. 

The  file  overflow,  in  (Appendix)  is  a  typical  input 
parameter  specification  of  the  cuci- symmetric  body  with  a 
conical  afterbody  using  the  grid  size  of  77x13x60.  This  file 
input  uses  NSTEP-100,  Mach  n\imber»0 . 95 ,  ALPHA«180°.  The 
calculation  method  depends  on  a  central  difference  Euler  term 
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in  J,  K,  L  and  ARC3D  diagonal  factorization  and  dissipation 
scheme  (IRHS-0,  ILHS»2,  IDISS*2) .  For  the  first  attempt  a  time 
step  DT-0.1,  ITIME=1  and  CFLMIN=5.0  was  utilized.  The  boundary 
condition  consists  of  6  boundaries  (NBC=6) ,  with  the  type  of 
BC  IBTYP=30  (inflow)  in  the  J  positive  (1),  IBTYP=15  (Axis  K 
round)  in  J  negative  (-1)  direction,  IBTYP=12  (symmetry  in  Y) 
in  the  K  positive  (2)  and  K  negative  (-2)  direction,  IBTYP=i 
(Inviscid  adiabatic  wall)  in  the  L  positive  (3)  direction  and 
IBTYP=32  (Supersonic  /subsonic  inf low/ out flow)  in  the  L 
negative  (-3)  direction. 

The  OVERFLOW  program  gives  output  files  such  as  ovr.out, 
q.save,  resid.out  and  fomo.out.  To  verify  that  a  calculation 
is  appropriate  or  converged,  it  can  be  traced  by  looking  at 
the  plot  of  residual  history  (resid.out) .  A  convergence 
criterion  was  defined  as  the  reduction  in  residuals  by  two 
orders  of  magnitude.  If  the  first  run  (NSTEP=100)  does  not 
fulfill  the  convergence  criterion,  one  may  perform  a  restart 
by  further  running  OVERFLOW  and  changing  the  parameter 
RESTART=.T.  in  the  input  file  (overflow. in)  and  by  copying 
q.save  into  q. restart.  Repeating  the  above  steps  until  all 
convergence  criteria  are  met.  If  in  the  output  file  q.bomb 
appears,  one  tries  another  run  by  changing  parameter  time  step 
DT  and  CFLMIN  until  the  criteria  are  met.  Then,  one  converts 
the  output  file  q.save  (unformatted)  into  formatted  file  named 
q.form  by  using  FORTRAN  code  READl.F  (Appendix) . 

Finally,  the  last  step  is  to  run  PL0T3D  program  with  the 
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source  formatted  files  x=grid.for  and  q=q.form.  By  using 
functions  on  the  PL0T3D  program,  many  plots  such  as  pressure, 
velocity,  Mach  number,  vorticity,  etc.,  are  obtained. 

3.  Results 

The  results  of  CFD  programs  are  grouped  into  the 
corresponding  geometric  shape,  namely,  afterbody  only  and 
complete  body.  The  Euler  flow- solver  (inviscid  flow)  was 
applied  to  all  axi- symmetric  bodies,  except  for  the  afterbody 
models  only,  where  both  Euler  and  Navier- Stokes  (viscous  flow) 
were  applied.  Most  of  the  calculations  converged  in  500 
iterations,  meaning  that  the  residual  history  achieved  a  two 
order  of  magnitude  drop. 

Results  can  be  analyzed  by  plotting  the  q.form  file, 
output  file  from  the  OVERFLOW  program,  using  the  proper 
function  in  PL0T3D  program  (Mach  number) .  By  interpreting  the 
Mach  number  contour  surrounding  the  body  surfaces,  one  can  be 
determine  the  characteristics  of  the  flow  field.  In  addition, 
the  drag  coefficient  C^j,  can  be  obtained  in  the  file  fomo.out, 
output  from  OVERFLOW.  Hence,  in  the  sequences  of  Mach  number, 
one  can  describe  the  significant  flow  characteristic  of  each 
geometric  shape. 

a.  Afterbodies 

An  Euler  as  well  as  a  Navier- Stokes  flow- solver  was 
applied  to  these  afterbody  models.  The  approaching  free  stream 
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Mach  number  (ranging  from  1.05  co  1.50)  starts  from  the  mid¬ 
section  of  the  afterbodies,  then  the  flow  follows  along  the 
afterbody  surface  until  it  reaches  a  maximum  local  Mach  number 
as  tabulated  in  Table  2 .  The  maximum  local  Mach  number  for 
viscous  flow  is  always  lower  compared  to  in  the  inviscid  flow. 
This  may  be  caused  by  the  viscous  flow  itself  since  we  are 
taking  into  account  the  shear  force  in  the  boundary  layer  near 
the  surface. 

For  the  afterbodies  from  the  small  perturbation  solution 
contour  (SPSl.O,  SPSl.l  and  SPS1.2) ,  in  inviscid  flow,  the 
shociis  are  formed  at  the  contour  surface.  The  location  of  the 
shock  depends  on  the  specific  afterbody  contour  and  the 
approaching  Mach  iiumber;  at  the  higher  Mach  numbers  the  shock 
appears  a  bit  further  downstream  as  shown  in  Figures  (16) 
through  (21)  .  As  shown  in  Figure  (22)  and  Figure  (23),  this  is 
apparently  a  result  of  viscosity;  it  shows  the  boundary  layer 
by  the  increment  of  Mach  nvimber  away  form  the  surface.  Weak 
shocks  were  formed  further  downstream  compared  to  inviscid 
flow.  Flow  separation  occurs  in  the  starting  contour  region 
and  is  followed  by  circulating  flow  in  the  base  region. 

For  the  conical  afterbody,  the  approaching  free  stream 
Mach  number  increases  following  the  mid- section  surface,  then 
a  flow  expansion  occurs  at  the  turning  angle  region,  until  the 
maximum  local  Mach  number  is  reached.  A  weak  shock  is  formed 
at  the  end  of  the  boat  tail  region  as  shown  in  Figure  (24)  , 
while  for  the  viscous  flow,  the  weak  shock  develops  away  from 
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boattail  surface  due  to  the  boundary  layer  and  the  separation 
of  the  flow.  The  circulating  flow  in  the  base  region  is  more 
significant  than  in  the  inviscid  flow  as  depicted  in  Figure 
(25) .  In  addition,  from  Table  2,  for  viscous  flow  the  maximum 
local  Mach  number  for  the  conical  afterbody  is  higher  than  for 
the  contoured  afterbody. 

For  each  afterbody  the  pressure  drag  coefficient  (C^j) 
versus  free  streeim  Mach  number  for  inviscid  and  viscous 
flow  are  plotted  in  Figure  (26)  .  The  negative  sign  of  C^j  is 
due  to  the  fact  that  the  calculation  of  pressure  starts  from 
mid- section  through  the  base  of  afterbody  and  ignores  the 
forebody  pressure.  These  results  show  that  the  pressure  drag 
coefficient  is  higher  for  viscous  flows  than  for  inviscid 
flows  for  each  given  afterbody.  This  may  be  caused  by  the 
viscosity  effect  and  the  pressure  distribution  difference  in 
the  flow  field.  In  addition,  it  can  be  seen  from  the  chart, 
that  the  C^j  for  the  afterbody  with  small  perturbation  solution 
contour  Mach  1.10  (SPSl.l)  has  the  lowest  values  over  the 
entire  Mach  number  range.  Therefore,  the  SPSl.l  contour  shows 
to  be  relatively  the  best  among  these  afterbodies. 

b.  Con^lete  Bodies 

Complete  bodies  consist  of  an  arbitrary  forebody 
(ellipsoid)  joined  to  various  afterbodies  such  as  the  small 
perturbation  solution  contours  (SPS_l.l  and  SPS_1.2)  and  a 
conical  afterbody.  The  approaching  free  streeun  Mach  number  (M^, 
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ranging  from  0.50  to  1.50)  starts  from  outer  boundary  with 
distance  of  1  times  for  the  coarse  grid  and  5  times  body's 
length  for  the  fine  grid  respectively.  Significant  differences 
between  coarse  and  fine  grid  calculations  can  be  seen  from  the 
characteristic  bow  shock.  In  the  coarse  grid,  the  bow  shock 
hits  the  outer  boundary,  this  causes  the  approaching  free 
stream  Mach  number  not  to  be  the  same  as  in  the  input 
parameter  as  shown  in  Figure  (27)  .  While  in  the  fine  grid,  the 
bow  shock  dies  out  before  reaches  the  outer  boundary  as 
depicted  in  Figure  (28) .  So,  further  discussion  is  focused  on 
the  fine  grid  exclusively. 

The  flow  stagnates  on  the  nose  and  then  follows  the  body 
surface  until  it  reaches  a  mcocimum  local  Mach  number  as 
tabulated  in  Table  3.  The  critical  Mach  nximber  for  these 
complete  bodies  is  approximately  at  M„  =  0.70,  where  the 
maximum  local  Mach  number  reaches  unity  at  the  shoulder 
region. 

Figure  (29)  shows  a  typical  high  stibsonic  free  stream  Mach 
number  (M^,  =  0.85)  flow  over  the  complete  body  SPS_1.1.  The 
flow  stagnates  on  the  nose  tip,  then  flow  is  accelerated 
following  the  forebody  surface  reaching  a  sonic  line  at  the 
mid-way  of  the  forebody  and  forming  a  supersonic  region  at  the 
shoulder.  Then,  the  flow  deccelerates  at  the  mid-section  and 
accelerates  again  until  reaches  a  maximum  local  Mach  number  at 
the  starting  region  of  the  afterbody  and  a  weak  shock  occurs 
in  this  region.  A  supersonic  approaching  free  stream  Mach 
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number,  is  depicted  in  Figures  (30)  and  (31)  for  SPS_1.1  cuid 
SPS_1.2  respectively.  A  bow  shocJc  is  obviously  seen  in  front 
of  the  nose.  The  characteristic  of  the  bow  shocJc  is  more 
inclined  down  stream  for  the  higher  Mo,  as  shown  in  Figure 
(32)  .  A  subsonic  region  is  formed  between  bow  shoclc  and  the 
nose,  then  the  flow  accelerates  along  the  forebody  surface  up 
to  a  supersonic  region  in  the  mid- section.  The  expansion  flow 
occurs  in  the  starting  contour  region  until  it  reaches  a 
maximum  local  Mach  niomber.  Then,  a  shoclc  is  formed  in  the 
contour  region.  Similar  as  in  the  afterbody  only,  the  shocJc 
location  depends  upon  the  contour  surface  and  Mj^.  The  shocJc 
location  for  a  given  contour  is  more  downstream  for  higher  M^ 
and  at  the  same  M„,  the  shocJc  location  for  SPS_l.l  is  more 
downstream  than  SPS_1.2.  A  typical  residual  history  of  CFD 
calculation  for  complete  body  with  small  perturbation  solution 
contour  is  shown  in  Figure  (33) .  Convergence  is  obtained  at 
about  500  iterations. 

For  the  conical  afterbody,  the  flow  characteristic  is  the 
same  as  the  other  complete  bodies  up  to  the  mid- section 
region.  The  expansion  flow  occurs  at  the  turning  angle,  then 
the  flow  accelerated  along  the  conical  surface  and  weaJc  shocJc 
is  formed  at  the  edge  of  base.  Figures  (34)  and  (35)  show  the 
Mach  number  contour  and  corresponding  residual  calculation  for 
conical  afterbody  at  M^,  =  1.10. 

The  pressure  drag  coefficient  (C^-press)  versus  free 
stream  Mach  number  (M*)  for  fine  grid  complete  bodies  are 
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plotted  in  Figure  (36)  .  The  drag  rises  sharply  in  the  high 
subsonic  Mach  number  (Mga  0.95)  and  reaches  a  meiximum  (pealc) 
at  «•  1.10.  Then  the  drag  decreases  with  a  shallow  curve  as 
the  Mj,  increases.  The  decreasing  shallow  curve  may  be  caused 
by  the  bluntness  of  the  nose  and  it  agrees  with  Shapiro 
[Ref .18]  because  the  fineness  ratio  and  bluntness  of  the  nose 
of  bodies  of  revolution  are  the  important  factors  that 
contribute  the  drag  curve  at  transonic  and  supersonic  range. 
It  can  be  seen  from  the  graph,  the  drag  curve  of  complete  body 
with  conical  afterbody  is  higher  than  with  small  perturbation 
solution  contour  at  the  entire  M^,.  Furthermore,  the  pealc  of 
the  drag  curve  is  approximately  15%  higher.  The  drag  curve  for 
SPS_1.1  and  SPS_1.2  are  lilcely  to  have  the  same  trend  up  to 
=  0.95;  beyond  this  Mach  number,  the  drag  curve  for  SPS_1.2  is 
slightly  greater  than  for  SPS_1.1.  Therefore,  the  complete 
body  with  small  perturbation  solution  contour  afterbody  Mach 
1.1  (SPS_l.l)  relatively  gives  the  lowest  drag. 
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VI.  CONCLUSIONS  AND  RSCOMMBNDATIONS 


The  use  of  numerical  simulation  (CFD)  appears  to  be  the 
most  cost  effective  method  to  predict  the  aerodynamic 
performance,  especially  in  the  transonic  range. 

In  this  research,  the  grid -generating  progreun  GRAPE  is 
suitable  only  for  the  geometry  with  smooth  contour  surface 
(SPS_1.1  and  SPS_1.2)  ,  while  the  prograun  GRIDGEN2D  is  used  for 
the  geometry  of  a  conical  afterbody  with  non- smooth  points. 

We  have  shown  that  for  a  complete  body  model  the  use  of  a 
fine  grid  (77x9x120)  is  more  reliadsle  than  a  coarse  grid 
(77x13x60)  .  This  is  shown  by  the  characteristic  of  the  bow 
shock  at  Mgs  >  1  as  seen  in  Figures  27  and  28.  The  pressure 
drag  coefficient  (C^j)  versus  free  stream  Mach  number  (M„) 
graphs  show  that  the  small  perturbation  solution  contour  for 
Mach  1.10  (SPS_1.1)  gives  relatively  the  lowest  on  both 
models  (afterbody  and  complete  body) ,  a  decrease  of 
approximately  15%  of  peak  in  the  transonic  range  compared  to 
the  conical  afterbody.  Therefore,  the  best  design  of  an  coci- 
symmetric  body  such  as  missiles,  projectiles  and  aircraft,  can 
be  based  on  the  small  perturbation  solution  contour. 

Finally,  this  work  may  be  continued  with  the  investigation 
a  complete  body  with  pointed  nose  geometry  and  the  application 
of  a  Euler  as  well  as  Navier- Stokes  flow- solvers  within 
various  small  angle  of  attack. 
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TABLE  1.  INPUT  PABAMBTBIt  FOR  OVERFLOW 
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TABLE  2.  MAXIMDIC  LOCAL  MACH  #  FOR  AFTERBODIES 


Inviscid  flow 


TABLE  3.  MAXIMUM  LOCAL  MACH  NUMBER  FOR  COMPLETE  BODY 
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HgiireS.  Numerical  integiatiaa  cf  equation  (63). 


SUPERSONIC  BOUNDARY  SURFACES 
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Figiire4.  Nvmetical  int^gratiaio  of  equation  (63)  (ii< 
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Figured  2D  grid,  afterbody  inoddia11iSPS_l.l 
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Hgnre?.  2Dgticl,  OQnq)letebod|fwit]iSPS_Ll. 
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Figures.  3D  grid,  afterbody  modd«attiSPS_l.l 
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Figure  10.  2D  grid,  conical  afterbody  model. 
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Figure  11.  2D  grid,  o(Mnplete  Ixxfy  with  conical  afterbody. 
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Figure  13.  3D  grid,  complete  body  with  ocmicalaftertxxly. 
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■grid,  oonq)Iete  body  with  SFS. 
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Figure  16.  Afterbody  model  SPS  l.l ,  invisdd  flow ,  Macb  1.10 . 
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Hgure  1&  Afterbod/ model  SPS_1.2 ,  invudd  flow ,  Macfa  LIO . 
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Figure  2L  Afterbody  model  SPS_1.0 ,  invisdd  flow ,  Madi  L20 . 
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Figure  23.  Afterbody  model  SPS  1.2 ,  viscouc  flow  ,  Mach  1 .20 . 
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Figure  24  Cboical  afterbody  model,  invisdd  flow,  Madi  1.10. 
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Figure  26.  ^  ^  ‘  Mach  number  for  aflerbocty  models. 
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Figure  28.  Fine-grid  (77x9x120)  for  quarter  complete  body  models. 
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Figure  29.  Quarter  complete  body  aft-SPS_l.l ,  Mach  0.85 . 
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Figure  30.  Quarter  complete  body  aft-SPS_l.l ,  Mach  1.10 . 
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Figure  31.  Quarter  complete  body  aft-SPS_1.2 ,  Mach  1.20 . 
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Figure  33.  Residual  history  of  complete  body  SPS_1.1 ,  Mach  1.10 . 
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Raaidual  Biatory 


ngure  35.  Residual  histoiy  of  ccMUi^te  body 


Cd. press  vs  Mach  *  for  complete  body  (fine  grid) 
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Figure  36.  -  press,  versus  Madi  number  for  quarter  complete  body  models . 


APPENDIX 

COMPUTER  PROGRAMS 


77 


PROGRAM  KSI.F 

•khli-itltitltlt***************************************************** 

*  This  program  is  developed  to  calculate  KSI (X) ,  K(X)  and  X  * 

*  using  mumerical  integration  based  on  trapezoidal  rule  to  * 

*  solve  equation  (32)  and  plotting  the  output  as  shown  in  * 

*  Figs .  1  .  * 

REALM(3),  X(0:401,3),  A,  Al,  B(401,3) 

#  ,  KSI,  P,  FUNG,  H,  DD,  Q,  L,  K(401,3),  N,  N1 

INTEGER  YY 

OPEN (UNIT=9 , FILE= ' KSI ' , STATUS= ' UNKNOWN' ) 

PRINT  *, 'ENTER  LOWER  &  UPPER  BOUND,  #  OF  INTERVALS,  #  OF 
#DATA  SET' 

READ  *,  A,  Al,  N,  N1 
100  DO  10  1=1,3 

PRINT*,  'ENTER  MACH  NO. ' 

READ*,  M(I) 

C  M(I)  =  0. 9+1*0. 1 

IF  (M(I) .LT.1.0)  THEN 
P  =  -1. 

ELSE 
P  =  1. 

ENDIF 

DO  20  J=1,N1 

B(J,I)  =  J*A1/N1 
H  =  (B(J,I)  -A)/N 
AREIA  =  0 . 

K(J,I)  =  (-*0.0102/(M(I)**2)  )*(  {ABS(  {B(J,I)**2)  -  (  (ABS{1- 
#  M(I)**2)  )  **1.7574)  ))**  (2. /3.)  +  (ABS(  1-M{I)**2)  )  **1.1716) 

c 

DO  30  L=  1,N 

KSI  =  A  +  (L-0.6)*H 

DD  =  (KSI**2+P*(ABS(1-M(I) **2) ) **1.7574) 

IF  (DD.GT.0.0)  THEN 
FUNG  =  1./ (DD** (1./3. ) ) 

Q  =  1. 

ELSE 

FUNG  =  l./( (-l.*DD)**(l./3.) ) 

Q  =  -1. 

ENDIF 

AREA  =  AREA  +  H*FUNC*Q 
X(J,I)  =  AREA 
30  CONTINUE 

20  CONTINUE 

10  CONTINUE 

C 

DO  40  J=  1,N1 

PRINT  50,  (B(J,I) ,X(J,I) ,K(J,I) ,1=1,3) 
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50  F0RMAT(1X,3  (F6.3,2F10.5,1X)  ) 

40  CONTINUE 

200  PRINT*, 'TRY  AGAIN?,  ENTER  1  FOR  YES,  2  FOR  DATA  FILE, 
#  OTHERS  FOR  NO' 

READ*,  YY 
IF  (YY.EQ.l)  THEN 
GO  'TO  100 
ELSE 

IF  (YY.EQ.2)  THEN 
DO  110  J-  1,N1 

WRITE  (9,50)  (B(J,I) ,X(J,I) ,K(J,I) ,1*1,3) 

110  CONTINUE 

GO  TO  200 
ENDIF 

ENDIF 

END 


PROGRAM  ZETA.F 

************************************************************ 

*  This  program  is  written  to  calculate  ZETA(r)  and  r  using* 

*  equation  (29)  and  plotting  output  in  Fig. 3  .  * 

********************************************* ♦***♦♦* ******** 

REALM(3),  ZETA(0;14,3) ,  R(14) 

OPEN (UNIT-9 , FILE- ' ZZ ' , STATUS- ' UNKNOWN' ) 

C 

DO  10  I-  1,3 

c  M(I)  -  0. 7+1*0. 1 

PRINT*,  'ENTER  MACH  NO.' 

READ*,  M(I) 

DO  20  J-  1,14 
ZETA(0,I)  -  1000.0 
PRINT*,  'ENTER  R  VALUES' 

READ*,  R(J) 
c  R(J)  =  0.1* J 

ZETA(J,I)  -  (4/(R(J)**2))+(ABS(l-(M(I)**2)) 

#  *(R(J) **2.8284) ) 

#  +( (l-M(I) **2)*(R( J) **7.657) /50. 63 

IF  (ZETA(J,I) .GE.ZETA(J-1,I)  ZETA( J, I) -ZETA( J-1, I) 

20  CONTINUE 

10  CONTINUE 
c 

PRINT  30,  M 

30  FORMAT (11X,5F10. 4) 

DO  40  J-  1,14 

WRITE(9,50)  (R(J),  (ZETA(J, I) , I-l, 4) 

50  FORMAT(4X,3 '2F10.4,2X) 

40  CONTINUE 
END 
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PROGRAM  D2D3.F 

*********************************************************** 

*  This  program  will  read  data  files  output  Gridgen2d  * 

*  from  2D  to  3D  for  further  processing  with  rotategr.f  * 

★★★★★★★★★A************************************************* 

REAL  X{77,60) ,Y(77,60) ,Z(77,60) 

READ (12,*)  IDIM,JDIM 

READ(12,*)  ( (X(I,J) ,I=1,IDIM) ,J=1,JDIM)  , 

#  ( (Y(I, J) ,I=1,IDIM) , J=l, JDIM) 

Z(I, J)=0.0 

WRITE (14,*)  IDIM,JDIM,1 

WRITE(14,*)  ( (X(I,J) ,I=1,IDIM) , J=l, JDIM) , 

#  ( (0.0,I=1,IDIM) , J=l, JDIM) , 

#  ( (Y(I, J) ,I=1,IDIM) , J=l, JDIM) 

STOP 

END 


PROGRAM  HALF.F 

*************************************************************** 

*  This  program  will  write  out  half  of  a  grid  file  output  * 

*  from  Grape  for  further  processing  with  rotategr.f  * 

It************************************************************ 

DIMENSION  X(200,200) ,Y(200,200) 

READ (30,*)  IDIM,JDIM 

READ (30,*)  ( (X(I,J) ,I=1,IDIM) , J=l, JDIM) , 

#  ( (Y(I, J) ,I*1,IDIM) , J=l, JDIM) 
c 

IDMD2=IDIM/2  +  1 
C 

DO  10  J=1,JDIM 

Y(IDIM2,J)  =0.0 
10  CONTINUE 
c 

WRITE (31,*)  IDMD2,JDIM,  1 

WRITE(31,*)  ( (X(I, J) ,I=IDMD2) , J=l, JDIM) , 

#  ( (0.0, I=IDMD2) ,J=1, JDIM) , 

#  ( (Y(I, J) ,I=IDMD2) , J=l, JDIM) 

STOP 

END 
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PROGRAM  ROTATEGR.F 

************************************************************* 

*  This  program  will  create  3D  grid  (body  of  revolution)  * 

*  by  rotating  3D  plane  grid  for  further  processing  with  * 

*  Flow  Solver  * 

************************************************************* 


DIMENSION  X (200, 200) ,  Y(200,200),  Z(200,200) 
DIMENSION  XX(200, 100,200) ,  YY (200 , 100 , 200) , 
#ZZ(200, 100,200) 

CHARACTER* 30  FNI 
CHARACTER* 30  FNO 

PRINT*,  'INPUT  GRID' 

READ (*,21)  FNI 

PRINT*,  'SG  =  0,  MG  =  1',  ',UNF-1,  FORM-2' 
READ(*,*)  IGRI 

PRINT*,  'IFORI  =  1  UNF,  IFORI  =  2  FORM' 

READ  ( * ,  * )  IFORI 
21  FORMAT  (A) 

IF  (IFORI. EQ.l)  THEN 
REWIND  1 

OPEN  ( 1 ,  FILE-FNI ,  FORM- '  UNFORMATTED '  ) 

IF  (IGRI. EQ.l)  READ  (1)  MGR 
READ  (1)  I1,J1,K1 

READ  (1)  ( (X(I,J) ,1-1,11) ,J-1,J1) , 

#  ( (Y(I,J) ,1-1,11) ,J-1,J1) , 

#  ( (Z(I,J) ,1-1,11) ,J-1,J1) 

CLOSE  (i) 

ENDIF 

IF  (IFORI. EQ. 2)  THEN 
REWIND  2 

OPEN  ( 2 ,  FILE-FNI ,  FORM- '  FORMATTED '  ) 

IF  (IGRI. EQ.l)  READ(2,*)  MGR 
READ(2,*)  I1,J1,K1 

READ (2,*)  ( (X(I, J) ,1=1,11) , J-1, Jl) , 

#  ( (Y(I,J) ,1-1,11) , J-1, Jl) , 

#  ( (Z(I,J) ,1-1,11) , J-1, Jl) 

CLOSE  (2) 

ENDIF 

PI  »  4.*ATAN(1.) 

PRINT*,  'NO  OF  PLANE  IN  J  DIRECTION  ?' 
READ(5,*)  JM 

DTH  =  (180./{JM-1) ) *(PI/180) 

DO  11  1=1,11 
DO  11  J-1,J1 
K=J 

XX(I,2,K)  =  X(I,K) 

YY(I,2,K)  =  0.0 
ZZ(I,2,K)  =  Z(I,K) 

11  CONTINUE 
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IM-Il 

KM-Jl 

DO  20  J-3,JM+1 
DO  20 

DO  20  K-1,KM 
XX(I,  J,K)-X(I,K) 

TH  -  (J-1)*DTH 
YY(I,J,K)  »  SIN(TH)*Z(I,K) 

ZZ(I,J,K)  =«  COS(TH)*Z(I,K) 

20  CONTINUE 

DO  30  I-l.IM 
DO  30  K«1,KM 
J-1 

XX{I,J,K)  =  XX(I,J+2,K) 

YY(I,J,K)  =-YY(I, J+2,K) 

ZZ(I,J,K)  =  ZZ(I,J+2,K) 

J-JM+2 

XX(I,J,K)  =  XX(I,J-2,K) 

YY(I,J,K)  --YY(I, J-2,K) 

ZZ(I,J,K)  =  ZZ(I,J-2,K) 

30  CONTINUE 
JM“JM+2 

PRINT'*,  'Output  filename 
READ (5, 21)  FNO 

REWIND  3 

OPEN (  3 ,  FILE= FNO ,  FORM= ' UNFORMATTED ' ) 

WRITE (3)  IM,JM,KM 

WRITE(3)  (((X(I,J,K),  I=1,IM) , J=l, JM) ,K=1,KM) , 

#  (((Y(I,J,K},  I=1,IM) , J=1,JM) ,K=1,KM) , 

#  (((Z(I,J,K),  I=1,IM) , J=l,  JM) ,K=1,KM) 
CLOSE  (3) 

STOP 

END 
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PROGRAM  READX.F 

It*********************************************************** 

*  This  progremi  will  read  3D  grid  file  and  convert  from  * 

*  formatted  into  unformatted  file  (grid. for  to  grid. in)  * 

*  for  further  processing  with  Overflow.  * 

*****************************************************it**i[ititi[ 

DIMENSION  X{ 77, 13,60) ,  Y(77,13,60),  Z(77,13,60) 

OPEN (UNIT- 12 , FILE- ' GRID . FOR' , STATUS- ' UNKNOWN' ) 
0PEN(UNIT-14, FILE- 'GRID -IN' , STATUS- ' NEW' , 

#FORM- '  UNFORMATTED ' ) 

READ  (12,*)  IDIM,  JDIM,KDIM 

READ(12,20)  (((X(I,J,K),  I-l, IDIM) , J-1, JDIM) , K-1, KDIM)  , 

#  (((Y(I,J,K),  I-l, IDIM) , J-1, JDIM) , K-1, KDIM) , 

#  (((Z(I,J,K),  I-l, IDIM) , J-1, JDIM) , K-1, KDIM) 

C 

WRITE (14)  IDIM, JDIM, KDIM 

WRITE(14)  (((X(I,J,K),  I-l, IDIM) , J-1, JDIM) , K-1, KDIM) , 

#  (((Y(I,J,K),  I-l, IDIM) , J-1, JDIM) , K-1, KDIM) , 

#  (((Z(I,J,K),  I-l, IDIM) , J-1, JDIM) , K-1, KDIM) 

20  FORMAT (6F15. 6) 

STOP 

END 


PROGRAM  READQ.F 

*********************************************************1,1,11 

*  This  program  will  read  q.save  file  output  from  * 

*  Overflow  and  convert  from  unformatted  into  formatted  * 

*  file  for  further  processing  with  Plot3d.  * 

************************************************************ 

DIMENSION  0(77,13,60,5) 

OPEN (UNIT-1 , FILE- ' Q . SAVE ' , STATUS- ' OLD ' , 

#FORM= '  UNFORMATTED '  ) 

OPEN  (20,  FILE- '  Q .  FORM'  ,  STATUS- '  NEW'  ,  FORM- '  FORMATTED '  ) 

READ(l)  NI,NJ,NK 

READ(l)  FSMACH,  ALPHA,  RE,  TIME 

READ(l)  ( (( (Q(I,J,K,NX)  ,I=1,NI) ,J=1,NJ) ,K=1,NK) ,NX-1,5) 
c 

WRITE (20,*)  NI,NJ,NK 

WRITE (20,*)  FSMACH,ALPHA,RE,TIME 

WRITE(20,*) ( ( ( (Q(I, J,K,NX) 


#,I=1,NI) , J-1,NJ) ,K=1,NK) , NX-1, 5) 
PRINT*,  NI,NJ,NK 
STOP 
END 
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PROGRAM  OVERFLOW.  IN 

•k*1fifk1t1t1t1t****1i*******1i1t1t***it***it****itiflc1t**-k*-k*1[**it*1i*1t*****1i 

*  Input  file  for  Overflow  to  run  either  EULER  or  NAVIER-  * 

*  STOKES  Flow- solver.  * 

************************************************************** 

$GLOBAL 

CHIMRA=.F.,  NSTEP=100,  RESTART=.F.,  NSAVE=100,  NQT=0, 
$END 
$FLOINP 

ALPHA=180.0,  FSMACH=0.95,  REY=  6.00E6,  TINF=  520.0, 
$END 
$GRDNAM 

NAME=  'AXI- SYMMETRIC  BODY  WITH  AFT- CONE,  77x13x60 
GRID' , 

$END 

$METPRM 

IRHS  =  0,  ILHS  =  2,  IDISS  =  2, 

$END 

$TIMACU 

DT  =  0.1,  ITIME  =  1,  TFOSO  =  1.00,  CFLMIN  =  5.00, 

$END 

$SMOACU 

ISPECJ  =  2,  DIS2J  =  2.00,  DIS4J  =  0.02, 

ISPECK  =  2,  DIS2K  =  2.00,  DIS4K  =  0.02, 

ISPECL  =  2,  DIS2L  =  2.00,  DIS4L  =  0.02, 

SMOO  =  1.00,  EPSE  =  0.35, 

$END 

$VISINP 

VISCJ  =  .F.  ,  VISCK  =  .F.  ,  VISCL  =  .F.  , 

NTURB  =  0, 

ITTYP  =  1 , 

ITDIR  =  3 , 

JTLS  =  1 , 

JTLE  =  77, 

KTLS  =  1 , 

KTLE  =  13 , 

LTLS  =  1 , 

LTLE  =  60, 

TLPAR1=  .3, 

$END 

$BCINP 


NBC 

=  6, 

IBTYP 

= 

15, 

12, 

12, 

1. 

32, 

15, 

IBDIR 

= 

1, 

2, 

-2, 

3, 

3, 

-1, 

JBCS 

= 

1, 

1, 

1, 

1, 

1, 

77, 

JBCE 

= 

1. 

77, 

77, 

77, 

77, 

77, 

KBCS 

= 

1, 

1, 

13, 

1, 

1, 

1, 

KBCE 

= 

13, 

1, 

13, 

13, 

13, 

13, 

LBCS 

s 

1, 

1, 

1, 

1. 

60, 

1, 

LBCE 

$END 

= 

60, 

60, 

60, 

1, 

60, 

60, 
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